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P. Argentiero 
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ABSTRACT 

In processing altimeter data from a spacecraft borne altim- 
eter to estimate the fine structure of the marine geoid, a 
problem is encountered. In order to describe the geoid fine 
structure, a large number of parameters must be employed 
and it is not possible to simultaneously estimate all of them. 
In practice, one is forced to hold a large number of param- 
eters at a priori values and adjust others. Unless the 
parameterization exhibits good orthogonality in the data, 
serious aliasing results. From simulation studies it has 
been found that amongst several competing par ameterizations, 
the mean free air gravity anomaly model (i. e. , Stokes’ 
formula) exhibited promising geoid recovery characterisitics. 

Using covariance analysis techniques, this paper provides 
quantitative measures of the orthogonality properties asso- 
ciated with the above mentioned parameterization. For in- 
stance, it has been determined that a 5° x 5° area mean free 
air gravity anomaly can be estimated with an uncertainty of 
1 mgal (40 cm undulation) provided that all free air gravity 
anomalies within a spherical radius of 10 arc degrees are 
simultaneously estimated. 
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STRATEGIES FOR ESTIMATING FROM ALTIMETER DATA 
THE MARINE GEOID 


INTRODUCTION 

The primary function of the spacecraft borne altimeter is to determine the fine 
structure of the mean sea surface topography. The instrument is well suited 
for the task. Consider, for instance, the altimeter which will be on board the 
GEOS-C spacecraft which is scheduled for launch in late 1974. In its global mode 
the instrument will be capable of producing a measurement every four kilometers 
along the subearth path of the satellite. This implies that even assuming con- 
siderable data compression it will be mathematically possible to extract from 
such data topographic detail of one degree or less. 

But there are practical difficulties to be overcome if the full potential of altim- 
etry as a data type can be realized. To see what these difficulties are it is nec- 
essary to closely analyze what happens when standard estimation techniques 
are employed to obtain sea surface topography from altimetry data. Essentially 
the problem is to reconstruct an analytic surface from direct but noisy observa- 
tions of the surface taken at specified points. The obvious approach is to param- 
eterize the surface in terms of coefficients which define a suitably dense set of 
functions in the space of two dimensional analytic functions and to recover the 
coefficients from the data by means of a standard minimum variance filter. Of- 
ten the physical setting of the problem suggests the proper parameterization. 

If not, then an arbitrary family of analytic functions such as the set of two di- 
mensional polynomials can be used. For this problem a natural parameterization 
is suggested by the fact that after suitable corrections the mean sea surface can 
be considered as cohering closely to an equipotential surface known as the marine 
geoid. (1) Thus one candidate for a natural parameterization is the set of stan- 
dard spherical harmonic coefficients of the Earth’s geopotential field. 

It is certain that any parameterization capable of describing the fine structure 
of the surface in question must employ an enormous number of coefficients. 

For instance, if we are interested in three degree features of the marine geoid 
and if the standard spherical harmonic expansion of the geopotential field is 
used as a parameterization, a full set of coefficients up to degree and order 60 
would be required. This implies the estimation of over 3700 parameters! Un- 
less very special circumstances apply it is not possible to simultaneously esti- 
mate such large parameter sets. In practice, in order to use parameter esti- 
mation techniques in recovering the marine geoid fine structure from altimeter 
data it will be necessary to adjust small subsets of the parameters while in effect 
'’freezing" the rest at a priori values. But unless the parameterization exhibits 
a certain property with regard to altimetry data which we term "orthogonality in 


1 



the data set", the net effect will be that the uncertainties in the adjusted terms 
will badly corrupt the estimates of the adjusted terms. This effect is frequently 
called aliasing. The orthogonality property just mentioned will be rigorously 
defined in a later section but in essence it is a property of a parameterization 
which permits a decomposition of the estimation problem into estimation prob- 
lems of much smaller dimensionality and without fear of serious aliasing. 

To take full advantage of the attractive properties of altimetry it will be necessary 
to develop a parameterization of the marine geoid which exhibits good ortho- 
gonality properties in altimeter data. The orthogonality properties of the set of 
spherical harmonic coefficients of the geopotential field is very poor and it is 
not a good candidate for the proper parameterization. Other candidates for the 
parameterization of the marine geoid are surface density layers (2) and sample 
functions (3). The parameterization whose properties will be investigated in this 
paper is the one provided by mean free air gravity anomalies and Stokes' formula 
(4), (5). If this parameterization is to be used to determine the marine geoid 
from altimeter data it will be important to discover to what extent mean free 
air gravity anomalies are orthogonal in the data. Specifically, we would like to 
know how far two mean free air gravity anomalies must be separated before one 
can be sure that a neglecting of one anomaly does not badly alias the estimate of 
the other. Without knowledge of this distance it is impossible to make intelligent 
use of the parameterization. 

The results of the SKYLAB altimeter experiment have demonstrated the ability 
of altimetry to reveal marine geoid fine structure. But it will be during the 
GEOS-C mission that for the first time scientists will have an opportunity to 
examine the output of a spacecraft borne altimeter functioning in a global mode 
over the world's oceans. In the succeeding section the performance of the 
SKYLAB altimeter and the goals of the GEOS-C experiment will be discussed. 
Next, the mathematical details of the measurement modeling and preprocessing 
of altimeter data will be described and the Stokes' formula mean free air gravity 
anomaly parameterization of the marine geoid will be documented. Following 
that will be a detailed discussion of the dual concepts of orthogonality and aliasing 
and their relationships to the problem of estimating a marine geoid from altim- 
eter data. Finally, the results of a systematic application of covariance analysis 
will be used to develop optimal estimation strategies for determining the marine 
geoid from altimetry data. 
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THE SKYLAB AND GEOS-C ALTIMETER EXPERIMENTS 

The GEOS-C spacecraft will be inserted into orbit during the latter part of cal- 
endar year 1974. This spacecraft will have as its prime experiment, a radar 
altimeter. The objectives of the GEOS-C altimeter experiment are: (a) to de- 
termine the feasibility and utility of a spaceborne altimeter to map the topography 
of the ocean surface with an absolute height accuracy of ±5 meters, and with a 
relative height accuracy of 1 to 2 meters, (b) to determine the feasibility of 
measuring the deflections of the vertical at sea, (c) determine the feasibility 
of measuring wave height, and (d) contribute to technology leading to a future 
operational satellite altimeter system with a 10 cm measurement capability. 

The altimeter data, when calibrated and corrected e. g. , for sea state, ocean 
tides, and other effects, constitute measures of the distance between the GEOS-C 
spacecraft and the ocean surface. Knowledge of the satellite altitude relative to 
a reference ellipsoid and knowledge of the oceanographic departures of the sea 
surface topography from the geoid will then permit the determination of the 
geoid. The chief problem is expected to be the determination of orbital altitude 
for GEOS-C. The primary tracking systems for doing this are the satellite to 
satellite tracking system and precision lasers. Data from these systems and 
others tracking GEOS-C, including C-band radars, USB range and range rate 
stations, and TRANET Doppler stations, will be used to find satellite height. 

Satellite contributions to the determination of the current ocean geoid have spa- 
tial resolutions corresponding to half wavelengths of approximately 18° (i. e. , 

2000 km). Surface gravity achieves representations with finer resolution, in 
the 1° to 5° range (i. e. , 110km to 550km), however it covers only about 50% 
of the ocean surface. The GEOS-C altimeter system will thus fill in the gaps 
and provide valuable independent knowledge where data now exist. 

The ability of a radar altimeter to detect features of the ocean surface has re- 
cently been demonstrated by the SKYLAB altimeter. In Figure 1 an altimeter 
pass over the Puerto Rican Trench has been analyzed. The pass was over the 
southwest corner of Puerto Rico and was 72. 8 seconds in duration. From a 
plot of the height residuals (formed by comparing the altimeter measurements 
with the calculated height of SKYLAB’ s orbit) it can be seen that a 17 meter 
drop in height residuals occurs when SKYLAB passed over the deepest part of 
the ocean (i. e, , corresponding to a 4000 meter depth of the ocean bottom) and 
the peaking of the height residuals occurred when SKYLAB traversed the Puerto 
Rican land mass. The SKYLAB data described here exemplifies the high level 
of resolution of surface features by a radar altimeter. 

It is to be pointed out however, that data from the GEOS-C experiment 
will be significantly free of the effects of spacecraft dynamic motions 
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which is not the case for SKYLAB altimeter data, since in the latter there were 
attitude control jet thrusting, crew motion, etc. In addition, because GEOS-C 
will have a nominal 1 year operating lifetime the altimetry data obtained will be 
far in excess of that obtained from SKYLAB. Thus for reasons cited above, the 
GEOS-C altimeter experiment data is most suitable for improving the marine 
geoid. 

To successfully utilize GEOS-C altimetry data for improving the accuracy of the 
marine geoid, dictated the development of unique computer programs capable of 
processing these data. The descriptions of the mathematical models associated 
with these programs and the computation strategies for processing altimeter 
data are provided in succeeding sections. 
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DESCRIPTION OF THE ALTIMETER MEASUREMENT 

The geometry associated with the altimeter measurement is described in Figure 
2. As can be seen, the altimeter measurement is nominally the shortest distance 
between the satellite and the sea surface, that is, the measurement is along the 
normal to the sea surface that passes through the satellite. The mathematical 
model for the altimeter measurement is given by the following relationship: 

h a = h-N'-5h s -Ah' (1) 


4 - 4 - , 

ip - r se l 

position vector of S/C 
geocentric radius vector 
geoid height above reference ellipsoid 
deviation of sea surface from geoid 

systematic errors in altimeter measurement (e. g. , refraction, 
timing, etc.) 

A more detailed description of the mathematical modelling of the altimeter meas- 
urement is given in Reference 6 (pp, 4-3 to 4-13). 

The error sources that effect the altimeter measurement fall into three cate- 
gories: those that are due to orbit uncertainty, those that cause the measured 
geoid to deviate from the true geoid, and those that effect the measurement ac- 
curacy itself. Equation 2 below describes the functional dependence of the 
error sources on the altimeter measurement: 

h a = h(E 0 t) - N'(G 0 ) - 5h s (? o ,l o ,0,\,t) - Ah r (m 0 ,^,X,t) (2) 


where 

h = 

4 

P = 


* 



Ah s 


where 

h: S/C altitude above reference ellipsoid 

E: orbit parameters and orbit dependent terms (radiation pressure, 

drag, etc.) 

t: time of altimeter measurement 

Cj 0 : vector of geopotential coefficients 
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Figure 2. Altimeter Measurement Geometry 









N f : geoid undulation function 

r 0 : vector of tidal error model coefficients 

s 0 : vector of local sea surface biases (i.e., currents, storm surges, 

wind waves. . . ) 

in 0 : vector of altimeter measurement error coefficients 
or in general 


— * “ > * 

h a ,Gq ,mQ (3) 

The altimeter measurement will be. corrected for known error sources, smoothed, 
sorted, etc. by the data preprocessor. Then, each measurement will be com- 
pared with altitude calculated from the satellite's orbit (h' •) to form the meas- 
urements residual Ah a : 


Ah a = h a -h' a 


Mathematically, this residual is equivalent to the difference function, 


3h • 9h a . 9h . 
Air = AE‘ + — AG J + — - Ar^ 


9F, 


1 

Bs 0 


9G; 


As k + 


J 

3ha 

Sim 




A m r + e 


or 


(4) 


(5) 


Ah a = Al^fE) + Ah a (G) + Ahgfr) + Ah a (s) + Ah a (m) (5a) 

The models for the altimeter measurement errors are now to be specified. 

1. Error Due to Orbital Parameters Ah a (E) 

Precision laser tracking data. Satellite to Satellite Tracking (SST), USB, C-Band 
and Doppler data preceding and during altimeter measurement data acquisition 
will be used to correct the error on altimeter measurements due to orbital pa- 
rameters using an existing orbit determination program. Orbit height accuracy 
required for altimetry data reduction must be better than 1 meter. Recent 
studies indicate that this is feasible (Ref. 7). 
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2. Error Due to Tides and Sea Surface Ah a (r) and Ah a (s) 


Corrections for deviations from the geoid due to ocean tides are based on the 
model of Hendershott (Ref. 8). In this model the dominant lunar semi-diurnal 
tide (M 2 ) with lunar declination terms neglected is represented. This tidal 
model is representative of the current state of the art. The maximum error 
contribution due tides is approximately 1 meter. The contributions due to local 
sea surface biases is on the order of 50 cm, below the resolution threshold of 
the GEOS-C altimeter. Therefore sea surface bias corrections are neglected. 

3. Errors Due to Altimeter Hardware Ah a (m) 

Altimeter measurement errors which directly affect the accuracy for determining 
the geoid radius are: 

Ah a (m) = Ah 0 +Ah£ + Ah L +Ah D +Ahj- +Ah R (6) 

where 

Ah Q : altimeter antenna offset 

Ahj: antenna offset due to S/C libration 

Ah L : dynamic lag 

Ah D : altimeter drift 

Ahj.: timing bias 

Ah R : tropospheric refraction 

Models for the above error sources are now described. 


A. Altimeter Antenna Offset 


Ah (J = h a tan Vt.{ol q - 5 b ) tan(a 0 - 6 b ) 


— 1 
.-—s 

N ' 
GO 

... 1 

Vi 

( a °\ - .225 

L 45 b h a_ 


_\ 25 b/ 


for 


5 b < «o < 2 ° 


(6a) 
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dr 



for 

0 < «o < 5 b 

where 

h a : spacecraft altitude 

5 b : antenna half beamwidth angle 

a 0 : angle off bote sight 

t : transmitted pulse length 

c: velocity of light (299.7925 x 10 3 km/ sec) 

B. Antenna Offset Due to S/C Libration 

Ahg = h a tan - 5 b ) tan(a e - 8 b ) 



CT) 

/ \ 

6, 


7 «o \ 2 . ' 

+ 



( — \ - .225 

1 

L 45 b h aJ 


Lw J 


for 

< «t < 2° 


or 


for 



0 < a £ < 5 b 


(6b) 


(6c) 


(6d) 
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where 


<*£ 

a c : S/C librati on angle 
: peak librati on angle 
<j>%: librati on phase angle 
T: orbital period. 

C. Dynamic Lag 

A measurement bias error induced by the change in S/ C position at the time of 
the outgoing and return pulses. 

9 *> 

Ah, = ! — - £ (6e) 

where 

o 

hji altitude rate at time t; (sec) 

X: lag coefficient (sec 2 ) 

D. Altimeter Drift 

A measurement bias error due to component aging. 

Ah D =D(t-t 0 ) (60 

where 

tg : initial time of altimeter measurement 
t: current time of altimeter measurement 

D: drift coefficient (sec) 

E. Timing Bias 

Bias error introduced by S/C clock error. 
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(6g) 


Ah T = hr 

o 

h: height rate at time of altimeter measurement 
r: clock error coefficient 

F. Tropospheric Refraction Error 

Bias error introduced into altimeter measurement due to tropospheric refraction. 
The Model used is that developed by J. Saastamoinen (Ref. 9). 

4. Error Due to Geopotential Ah a (G) 

In the preprocessing of altimeter data the altimeter measurement will have been 
corrected for the orbit uncertainties ocean tides, local sea surface biases and 
measurement bias errors using the error models described above. Thus Equa- 
tion 4 above is restated as follows: 

Ah^G) = h a - [h; + Ah a (E) + Ah a (r) + Ah a (s) + Ah a (m)] (7) 

And the relationship between altimeter measurement residuals to geoidal param- 
eters can then be stated as follows: 

9h a 

Ah a (G > “-£7 & G j - < 8 > 

SCj 

As stated above, the altimeter measurement residual as described in (7) solely 
reflects primarily (neglecting second order effects) the departure of the geoid 
from the reference ellipsoid. Or more precisely the distance along the UNIT 
normal to the reference ellipsoid between the reference ellipsoid and the geoid 
which is called geoidal height or geoidal undulation N' (see Figure 2). 

The mathematical relationship which relates the geoidal undulation to the dis- 
turbing potential T is given by Bruns' formula derived in Reference 4. That is 

T 

N =— (9) 

y 


where 

T = disturbing potential 

7 = magnitude of gravity vector normal to reference ellipsoid 
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The disturbing potential of the global geoid at point (<j>, X, r) is expressed in terms 
of spherical harmonics as follows: 

L „ n 
GM t — ^ /a\ ▼ — * 

T S (7 ) 2 P .m< sin f C nm cos mX + S nm sin mX J < 10 > 

r n= 2 V ' m = 0 

where 

GM: gravitational constant of the earth 

C nm , S nm : spherical harmonic expansion coefficients 

P |im (sin 4>)\ associated legendre functions 

(0,X): latitude and longitude on geoid at which disturbing potential is 

evaluated 

r: geocentric radius from earth center of mass to evaluation 

point of disturbing potential on geoid 

a: semi-major axis of reference ellipsoid 

L: is the limit of summation, and it is specified by the degree of 

harmonic expansion of the global geoid 

n: summation index for degree terms of the spherical harmonic 

expansion of T 

m: summation index for the order of terms in spherical harmonic 

expansion of T 

Thus the geoidal undulation at any point P (<p,\, r) on the earth can be computed 
from geopotential coefficients derived from satellites by analysis of perturbations 
on the orbits induced by the earth's gravity field. The undulations are computed 
from the combination of Equations 9 and 10. 

Another form of expressing the disturbing potential is in terms of Stokes' formula 
(Ref. 4 pp. 92-98). This formula makes it possible to express the disturbing 
potential of the global geoid from gravity data. That is 

T = £ J f AgS(^)d a (11) 

•'a * 
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where 



R = a( 1 - 0 

a: 

radius of reference ellipsoid 

f: 

flattening of reference ellipsoid 

S(0) = 

Stokes' function 

a: 

element of area 

Ag: 

surface gravity anomalies 


S<0) = esc 


/ \b \ i jj 

\2/ ~ 6 ?in 2 + 1 - 6 cos ^ " 3 cos 




(11a) 


from Bruns' formula (Eq. 9). 

The geoidal undulation at any point P on the global geoid can be computed from 
Stokes' formula (Eq. 11). That is 



R 

47T7 



Ag S(i //) do 


(Hb) 


In terms of geographical coordinates Stokes* function can be expressed as 
follows: 


where 


N(0,X) = 


R 

47T7 


27 r r 7r/2 


J Z7T /*7 


Ag(0',X') S(0) cos b'd0'dX' 


(He) 


X' = 0 V =■-»/ 2 


do = cos0'd0'dX' 

(0,X) = latitude and longitude of the computation point 

(0'X') = coordinates of the variable surface element a 

0: spherical distance between the computation point and variable 

surface element 


0 = cos 1 


[sin 0 sin <p' + cos <j> cos <j>' cos (X - X r )J 
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Ag(0'V): free air gravity anomaly at the variable point (<p'X') 

Ag (0'V): fe ( rtF riable pomt ( ^ X * 

y- mean value of gravity over earth ...» . . , . , 

In order to combine surface gravity data and geopotential information derived 

the point JP. Each gravity ^om^.m each area ^g c pfetbpftrtati£m^4ht0^®)pai^ 



represented by _A^ and^Ag 2 respectively where 


Ag. = 7 


Ag s = 7 


(lid) 

(lid) 


"£ X (n " 1 } (“)n P nm (sin & { C nm cos nX + S nm sin mX > , 
n^} n^O (n - 1) t'A P nm (sin d)f nm cos nX +S nm sin mX}J 

The A g 2 value i^ 1 defied 2s the remainder of the gravity anomaly. By parti- 

the v |^h-^g f ^ d a£^^g^^^ r g^c^ltM<^tetiai?sahd3Qn.!tl?Q c&yLpb^irts. 



»mng 


value is aennea se tuct ^ ^ ^ ~ * 

!te nl undulations into two components 

Eauation 11c can thus be rewritten as follows: 

N(0,X) = N 1 +N 2 (lie) 

N(A,X) = N,+N 2 0 1cX 


where 

where 


-2t r »ir/2 


N, = 


or 

or 

and 

and 


N ‘ = L2 tt| »'/r/lAg s (0',X')S(i//) cos ^'d^'dX' 
jvj = '-|r/2 Au.(^ , ,X , )S(^)cos^'d0'dX' 

1 4 ^' T 0 K n 


Z 

7r 

y 


R 

N» = 

2 47Tt- 



n 2 = ~r \ 

2 4*7 J 


/Ag 9 (0',\') S(^) cos 0 r db'dX' 

* Ag, ib',A')S(- v ■ osd'dd'dX 


A J 


~ AT 

N 2 = ANj = 5N^= 

N 2 = AN, = 6N =TA- 

Ag 2 = 5(Ag s ) 

Ag-, = §(Ag : 
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Thus 


N(0,X) 


T + AT 
7 


n 2 = an l 


4tt7 

j=l 


5{Ag s (0]X])}S(V/j) cos 0- 


(Ilf) 


where ANj is the correction to the geoidal undulations of the global geoid as a 
function of. the corrections of mean free air gravity anomalies. 


Equation Ilf is the form of the- parameterization adapted for relating the altim- 
eter measurement residuals to geoidal parameters. That is. Equation 8 can 
now be written as fellows: 


5N = Ah a = ANj 


or in the matrix form 


§g = AGj = 5{Ag s (0- Xj)} 

3h, R A<b'AX' 

A = = S(i^.) cos-^j 

0G| 4^7 J J 


^ N (k X 1) A (k X j) 5 S(jX l) 


( 12 ) 


(12a) 
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ORTHOGONALITY AND ALIASING 


Assuming that altimeter observations can be directly translated into deviations 
of the marine geoid from a reference geoid, it^s straightforward to estimate 
gravity anomalies from altimetry data. Let 5N be a vector of geoid undulations 
collected from a certain region over the ocean. Next let 5g be a vector of mean 
free air gravity anomalies defined over a region of the ocean which contains the 
data region. Then repeating Equation 12 we have 

5N = A 5g (13) 


where A is a matrix the number of whose rows is equal to the number of data 
points and the number of whose columns is equal to the number of gravity anom- 
alies. As demonstrated in Eq. 12 the individual elements of A, say A (I, J), can 
be computed through Stokes' formula and the latitude and longitude of the i data 
point together with the longitude and latitude of the midpoint of the grid over 
which the J th gravity anomaly is defined. Equation 13 provides a linear equation 
of condition and in a standard minimum variance fashion 5g could be estimated 
from observations of 5N. But in order for Equation 13 to be correct (correct, 
that is, assuming that the approximations inherent in the discrete form of Stokes' 
formula are valid) the gravity anomalies in the array Sg must cover the globe. 
Considerably smaller regions would no doubt be adequate but just how small 
these regions can be before a serious bias is introduced into the estimation 
process is a matter to be determined from computations. In any case, compu- 
tational consideration constrain us to chose a region such that the number of 
elements in the Sg array does not exceed two or three hundred. Gravity anom- 
alies outside of this region are in effect assumed to be zero. To see precisely 
what happens when this assumption is made, postulate that the Sg array of 
Equation 13 is defined over an area sufficiently large that Equation 13 is sub- 
stantially correct. Then write 


8g = 


5g 2 


(14) 


where 

5gj = gravity anomalies to be adjusted in a standard minimum variance 
filter 


and 

6g 2 = gravity anomalies assumed to be zero and thus left unadjusted by the 
minimum variance filter. 
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Then Equation 13 can be written 

5N = A 1 Sgj + A 2 5g 2 (15) 

where A, and A 2 are respectively the variational matricies of 8N with respect 
to 8gj and with respect to §g 2 . 

After proper corrections altimeter data is treated as direct observations 5N of 
5N with statistics 

5N = 5N + v, E(r) = 0, H(w T ) = Q (16) 

Estimates of mean free air gravity anomalies obtained from satellite tracking 
and gravimetry measurements are available. 10 Unless this information were 
correctly factored into the gravity anomaly estimates obtained from altimeter 
data, the resultant estimates would not be optimal. Consequently- we assume 
the existence of an a priori estimate 5gj of 5gj with properties 

5 g ; = 5g, +<*,,£(«,) = 0,E(a,«T) - P, (17) 

The gravity anomalies Sg 7 for computational reasons are assumed to be zero 
but the actual values of gravity anomalies in the region on the sphere which is 
ignored have a certain distribution about zero. We assume 

E(Sg 2 ) = 0, E(5g 2 6g 2 T ) = P 2 (18) 

When the values of 8g 2 are assumed to be zero the minimum variance estimate 
of 5gj becomes 


8gj = (A} Q- 1 A, +P- 1 )- 1 [A} Q' 1 5N + P- 1 5 g ; ] (19) 

See Reference 11 for details. Define the covariance matrix of the estimator 
given by equation 19 as 


P = E [(Sg, -SgjHSgj -§ gl ) T ] (20) 

From Equations 15, 16, 17 and 19 we obtain 

5gj -5g L = (A{ Q" 1 A t +P7 1 )" 1 (~A| Q- 1 A 2 8g 2 + A} Q~ l V+?" 1 a { ) (21) 

Equation 21 yields 


P = (A| Q -1 Aj +P 7 1 )- 1 +(A[Q- 1 A,Pj 1 r 1 A| Q" 1 A 2 P 2 A 2 Q _1 A, 


(A{Q' 1 A 1 + P] 1 ) 


( 22 ) 


lv-i 
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Assume that the data is uncorrelated and that each data point has the same 
variance. Then 


Q = I*o (23) 

Where I is the identity matrix and a 2 is the common variance of the data. Also 
assume that the values of the unadjusted gravity anomalies are independently 
distributed. Then the covariance matrix P 2 of 6 g 2 can be written as 



Where n 2 is the number of unadjusted gravity anomalies and a 2 is the second 
moment about zero of the i 1 * 1 unadjusted gravity anomaly. • 

Also define a matrix K as 

K s (AjQ^Aj +P 7 1 )- 1 A}Q- 1 A 2 (25) 

If n t is the number of adjusted parameters, then K is of dimension by n 2 . 
With these assumptions Equation 22 yields the following expression for the 
variance of the i tfl adjusted gravity anomaly 

n 2 

P(U) - £ W W] up 2 (26) 

j = 0 

Where |0 2 o is the i th diagonal element of the matrix (A]A 2 ) -1 (we assume here 
that diagonal elements of the matrix Pj 1 are relatively small) and 

0ij = K(I,J), J > 1 (27) 

Define the error sensitivity matrix as 

s — > ( ~ J — 0,1, 2,.,,. n 2 (28) 

and finally define the Alias Matrix as 


L = so 


(29) 



where 



0 



( 30 ) 


The alias matrix reveals much of the probability structure of the estimation 
procedure. From Equations 26, 28, 29 and 30 it can be seen that the standard 
deviation of the i th adjusted gravity anomaly is the root sum square of the terms 
in the i 1 row of the alias matrix. The elements in the first column of the alias 
matrix represent the R. S. S. contribution to the standard deviation of each esti- 
mated parameter due to the data noise. The elements in the J tln column, J > 2, 
represent the R. S. S. contribution >to the standard deviation of each estimated 
parameter due to the J - 1 st unadjusted parameter. These terms are called 
the aliasing contributions to the uncertainty in the adjusted parameters due to 
the uncertainty in the J - 1 st unadjusted parameter. Notice that the aliasing 
contributions due to the J ^ unadjusted parameter are proportional to the stan- 
dard deviation of the J* parameter. Definition: hiagiven estimation process the 
i 1 estimated parameter is said to be orthogonal with respect to the J th unesti- 
mated parameter if the aliasing contribution to the i * ^ estimated parameter due 
to the uncertainty of the J 1 ' 1 unestimated parameter is zero. 

Two things can be noticed concerning this definition. First, the orthogonality 
relationship between two parameters must be defined within the context of a 
given estimation process. The data set must be defined and it must be stipulated 
which parameters are in the adjusted and in the unadjusted mode. Second, al- 
though we have been discussing mean free air gravity anomalies and altimeter 
data our mathematical results and definitions are applicable to any linear esti- 
mation problem where some parameters are adjusted and others left unadjusted. 

To see the implications of the orthogonality relationship we need a more re- 
vealing representation of the aliasing terms. First notice that the first term 
on the right side of Equation 22 is the covariance matrix of the estimation 
process under the assumption that the unadjusted parameters are perfectly 
known. This covariance matrix gives the uncertainty of the estimates due only 
to data noise. Define the so-called "noise only" covariance matrix as 

P s (A[Q- 1 A 1 +P7 1 )- 1 (31) 
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Next observe that the elements in the i 1 * 1 row and J 1 * 1 column of respectively 
Aj and A 2 are the partial derivative of the i t(l data point with respect to the 
J th adjusted parameter and the partial derivative of the i th data point with re- 
spect to the J tfl unadjusted parameter. The aliasing contribution to the i lil ad- 
justed parameter due to the J th unadjusted parameter can be written as 


n t m 

m,j + 1) = £ Pd,K) 

K = 1 L = 1 

™ Q-.(L,L) ™ 
38 gl (K) ’ 36g 2 (J) 

where m is the number of data points. 


(32) 


If the estimates of the adjusted parameters are relatively uncorrelated in the 
noise only covariance matrix. Equation 32 can be approximated by 


m,j + 1) - p(i,D 


V d g N(L) t 85N(L) 
35^(1) g (L ’ L) 85g 2 (J) 


(33) 


A sufficient condition for the left side of Equation 21 to approximate zero is for 
the observability patterns of 5g 1 (I) and Sg 2 (J) in the data to be virtually non- 
overlapping. Up until about 40°, Stokes’ function rapidly attenuates with in- 
creasing values of spherical radius. 4 Hence if the grids on which 6g 1 (I) and 
5g 2 (J) are defined are sufficiently separated, the orthogonality relationship 
would be effectively satisfied and the estimate of Sg ^ (I) would experience no 
aliasing from the uncertainty of 5g 2 (J). Conversely if the grid on which 5g 1 (I) 
was defined were in close proximity to grids whose gravity anomalies were un- 
adjusted, one would expect serious aliasing of the resultant estimate. 

It should be clear then, that if the gravity anomalies are estimated in a block 
the outer layers of the block contain gravity anomalies whose estimates will be 
badly aliased by the adjacent unadjusted parameters. It will be necessary to 
discard these estimates. But the gravity anomalies in a sufficiently small inner 
core of the block may be adequately separated from the unadjusted parameters 
to be effectively orthogonal with respect to them. The estimates of these terms 
presumably will be of sufficient accuracy that they can be accepted. In effect, 
for every block of gravity anomalies that we intend to estimate it will be neces- 
sary to construct a ’buffer zone” several layers deep of gravity anomalies which 
surround the block. The new and larger block of gravity anomalies must be 
simultaneously estimated and then the estimates of gravity anomalies in the 
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buffer zone must be rejected due to aliasing. In order to achieve an intelligent 
compromise between estimation accuracy and computational load it is necessary 
to know the relationship between the depth of the buffer zone and the accuracy of 
the estimation procedure. The relationship will vary with grid size, computation 
radius and data density. The most convenient way to study the relationship is 
to utilize covariance analysis techniques to generate alias matricies for several 
situations and to attempt generalizations from the results. This is accomplished 
in the following section. 


22 



ESTIMATION STRATEGIES FOR GRAVITY ANOMALY DETERMINATION 


The accuracy with which a given gravity anomaly is estimated from altimeter 
data is a function of many parameters. It is of course dependent on the accuracy 
and density of altimeter data. As explained in the previous section it is also 
dependent on the radial distance between the estimated gravity anomaly and the 
nearest gravity anomaly which is in the unadjusted mode. This parameter we 
call the estimation radius. Another parameter, the computation radius, is an 
important determinant of the accuracy of a gravity anomaly estimation. The 
computation radius sets the maximum distance from a given grid element over 
which data is to be processed in order to estimate the gravity anomaly defined 
on that grid element. Estimation accuracy does not necessarily increase with 
increasing computation radius. To see the reason notice that for a given esti- 
mation radius the covariance matrix of a set of estimated gravity anomalies is 
given by Equation 22 as the sum of a matrix which is dependent only on data 
uncertainty and a matrix which represents the aliasing effects from the unadjusted 
parameters. With increasing computation radius the elements of the first matrix 
must decrease but in general the elements of the matrix which conveys the ali- 
asing effects will increase. This effect can be shown graphically by means of so- 
called aliasing maps. To obtain our aliasing maps of Figures 3 and 4, we simu- 
late in a covariance mode a situation in which we describe the marine geoid by 
means of two degree by two degree gravity anomalies and we assume 12 alitmeter 
observations for each grid, each with an uncertainty of 1 meter. In Figure 3 
we map the R. S. S. contribution in mgals to the uncertainty in the estimated 
gravity anomaly defined on the grid element in the lower left hand corner when 
the adjacent anomalies are assumed to be in an unadjusted mode and to have an 
a priori uncertainty of 50 mgals. The computation radius is 5°. Notice that 
the aliasing contributions decrease with the distance between the unadjusted 
parameter and the estimated parameter thus demonstrating the inherent ortho- 
gonality property of the gravity anomaly parameterization of the marine geoid. 

In Figure 4 the computation radius has been changed to 20° and the aliasing effect 
is seen to decrease much less radically. 

To achieve an intelligent compromise between computational load and estimation 
accuracy it is necessary to determine the precise relationship between the ac- 
curacy of the estimate of a given gravity anomaly and the estimation radius and 
computation radius employed in the estimation procedure. To do this we assume 
that altimeter data existed at a density of three data points per one degree by 
one degree grid and that the uncertainty on the data was one meter. It was also 
assumed that unadjusted gravity anomalies had a standard error about zero of 
50mgal. This figure is conservative, a more realistic figure being 30mgal 12 . 

No a priori estimate was assumed for estimated parameters. The effect of 
unadjusted parameters out to a spherical radius of 45° were included in the com- 
putations. Under these conditions Figure 5 provides the standard deviation in 
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Figure 3. Alias Map for 2° by 2° Grids 
and for 5° Computation Radius 
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Figure 4. Alias Map in mgals for 2° by 2° Grids 
and for 20° Computation Radius 
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COMPUTATION RADIUS 


Figure 5. Accuracy of 5° by 5° Mean Free Air Gravity Anomaly Estimate vs. Computation Radius 

for Various Estimation Radii 



the estimates of 5° by 5° mean free air gravity anomalies as a function of esti- 
mation radius and computation radius. From this figure it is clear that the most 
efficient estimation strategy for 5° by 5° anomalies is obtained with an estimation 
radius of about 10° and a computation radius of about 5°. This strategy will re- 
sult in an accurac}' of 1 mgal. This implies that if 5° by 5° mean free air gravity 
anomalies are estimated in blocks, estimates of gravity anomalies in the outer 
two layers of the block should be discarded due to aliasing. The rest of the 
estimated gravity anomalies should be accurate to within one mgal provided a 5° 
computation radius is used. It is important to remember, however, that our 
prime intent is to estimate a marine geoid rather than gravity anomalies. With 
estimates of 5° by 5° mean free air gravity anomalies can utilize Stokes' formula 
to analytically reconstruct a marine geoid which is sufficiently detailed that 5° 
features can be noticed. Since the discrete form of Stokes' formula is linear it is 
an easy matter to propagate a given error in gravity anomaly estimation into an 
error in marine geoid determination. Assuming that Stokes' formula is accurate 
if its summation is carried out within a 90° spherical radius of a given point, a 
one mgal standard deviation in the estimate of 5° by 5° mean free air gravity 
anomalies propagates into a 40 cm standard deviation of the resultant marine 
geoid. Thus by applying proper estimation procedures to the reduction of altim- 
eter data should be possible to determine 5° features of the marine geoid with a 
resolution of 40 cm. 

To obtain a more detailed marine geoid represents a more difficult estimation 
problem. Figure 6 provides the standard deviation in the estimates of 3° by 3° 
mean free air gravity anomalies as a function of estimation radius and computa- 
tion radius. An estimation radius of 12° and a computation radius between 10° 
and 11° appears to provide the most intelligent estimation strategy. The strategy 
leads to estimates of 3° by 3° mean free air gravity anomalies which have stand- 
ard deviations of 5 mgal. This implies that is 3° by 3° mean free air gravity 
anomalies are estimated in blocks, the outer 4 layers must be discarded as 
being badly aliased. A 5 mgal standard deviation in estimates of 3° by 3° mean 
free air gravity anomalies propagates into a standard deviation of 1.2 m in the 
resultant marine geoid. With the present assumptions no estimation strategy 
is adequate to determine features of the marine geoid as small as 2°. With much 
greater data densities such fine resolution is possible. For instance if we increase 
the data density by a factor of 100 thus postulating 300 data points for each 1° by 
1° block, then with an estimation radius of 10° and a computation radius of 2°, 

2° by 2° anomalies can be estimated with a standard deviation of 2, 3mgals. This 
propagates into a standard deviation of the marine geoid of about 42 cm. But to 
postulate such data densities is probably not realistic and in addition it suggests 
a very heavy computational load. The accurate resolution of 3° features of the 
marine geoid is likely to be the limit of what can be accomplished with altimeter 
data and the Stokes' formula mean free air gravity anomaly parameterization of 
the marine geoid. 
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Figure 6. Accuracy of 3° by 3° Mean Free Air Gravity Anomaly Estimate vs* Computation Radius 

for Various Estimation Radii 




Finally notice that these results were accomplished without assuming a priori 
estimates of gravity anomalies. In fact such estimates have been derived from 
satellite tracking data and direct gravity measurements. When these estimates 
are optimally combined with altimeter data by means of Equation 10 the resultant 
estimates of gravity anomalies will be predicted on all relevent data. The repre- 
sentation of the marine geoid derived from these estimates will then constitute 
the best possible estimate of that surface. 
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SUMMATION 


Data from the SKYLAB altimeter has demonstrated the ability of altimetry to 
disceme the fine structure of the marine geoid. It will be during the GEOS-C 
mission, however, that global sensing of the sea surface topography from a 
spacecraft borne altimeter will be possible for the first time. In order to make 
optimal use of the vast amomit of altimeter data expected from the GEOS-C one 
must apply numerous corrections to the data, both to eliminate systematic errors 
inherent in the data type itself and to correct for the effects of deviations of 
the mean sea surface from the marine geoid. 

In order to employ standard parameter estimation techniques to estimate a 
marine geoid from altimeter data it is necessary to utilize a parameterization 
of the marine geoid which exhibits good orthogonality characteristics in the data 
type. In this respect at least, the Stokes' formula mean free air gravity anomaly 
parameterization appears to be the most promising. For a given grid size and 
data density, the accuracy of the geoid resulting from an estimate of gravity 
anomalies from altimeter data is a strong function of the choice of estimation 
radius and computation radius. An application of covariance analysis techniques 
reveals that assuming a data density of three measurements per one degree 
by one degree spherical surface area with each measurement accurate to within 
one meter, optimal choices of estimation radius and computation radius of 
respectively 10° and 5° yield an accuracy in the estimate of 5° by 5° mean free 
air gravity anomalies of approximately 1 mgal. This result implies that with a 
judicious choice of estimation strategy, one can determine geoid height with a 
standard deviation of 40 cm and a spatial resolution of 500 km (5 arc degree). 

With the same data density and data accuracy assumptions, covariance analysis 
demonstrates that optimal choices of 12° for estimation radius and 10 for 
computation radius permits one to estimate with an accuracy of approximately 
1. 2 m in geoid height with a spatial resolution of 300 km (3 arc degrees). Severe 
aliasing difficulties are encountered in attempting to estimate more detailed 
geoids. Our covariance analysis studies indicate, however, that a very small 
computation radius together with great data densities can mitigate the aliasing 
effects and yield in localized areas a marine geoid capable of showing spatial 
resolutions at the 100 km to 200 km level. 

Finally it should be mentioned that by utilizing the Stokes' formula mean free 
air gravity anomaly parameterization of the marine geoid, it becomes an easy 
matter to combine satellite tracking data and direct measurements of gravity 
anomalies with altimeter data to obtain an optimal estimate of the marine geoid. 
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